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Abstract 

'a t 

We model the density dependence of the Griineisen parameter as j(p) = 1/2 + 
^ . Ti/p 1 ^ 3 + 72/ P 9 , where 71, 72, and q > 1 are constants. This form is based on the 

assumption that 7 is an analytic function of V 1 / 3 , and was designed to accurately 
represent the experimentally determined low-pressure behavior of 7. The numerical 
values of the constants are obtained for 20 elemental solids. Using the Lindemann 
criterion with our model for 7, we calculate the melting curves for Al, Ar, Ni, Pd, 
^ ■ and Pt and compare them to available experimental melt data. We also determine 

the Z (atomic number) dependence of 71. The high-compression limit of the model 
is shown to follow from a generalization of the Slater, Dugdale-MacDonald, and 
Vashchenko-Zubarev forms for the dependence of the Griineisen parameter. 
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Introduction 

Lattice anharmonicity leads to a volume dependence of the phonon frequencies, Ui, that 
is described by the mode Griineisen parameters JTfl 

ain^ 

ll dlnV 1 j 
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and the mechanical Griineisen parameter is defined as the average of the 7, over the first 
Brillouin zone, 7 m = (7$). When the % are all equal it can be shown that 7 m = 7$ 
coincides with the thermodynamic Griineisen parameter |J 

1th = — , (2) 

where a is the thermal expansion coefficient, Bt is the isothermal bulk modulus, and Cy 
is the heat capacity at constant volume. If the 7, are not all equal, 7 m 7^ jth m general. 
As noted by Barron ||, Born has shown that, even if the ji are not all equal, 7 m = 7^ in 
the limit of low (T — - > 0) and high (T ~ 0£>) temperatures. The T = limiting value of 
7m = 7t/i = 7 is 7 = — dln©£>/G?ln V, where G>£> is the Debye characteristic temperature, 
^(u;^ 3 ) -1 ' 3 /^, in the limit T -> |. Here both 0£> and 7 are functions of volume, 
or density, alone. This formula is the Debye-Griineisen definition of 7 [|J, which can be 
rewritten in terms of density, p ~ 1/V, as 

rfln6 D (p) 

7(P) = — dl^— " (3) 

The Lindemann melting criterion, which asserts that the root-mean-square atomic 
displacement of atoms from their equilibrium positions in a solid is a fixed fraction of the 
interatomic distance at the melting point, can be rewritten in the form of the Gilvarry law, 
which relates the density derivative of the melting temperature, T m (p), to the Griineisen 
parameter M: 



d\nT m (p) 



7(P) " ~ 



(4) 



dlnp 

hence integration of 7(p) yields the melt curve. Equations of state for solids can be 
constructed on the basis of y(p), though additional thermodynamic data are needed ||. 

There is a long history of attempts to model y(p) J7|. Slater [|J, Dugdale and Mac- 
Donald ||, and Vashchenko and Zubarev |T(| proposed three expressions for j(p) that 
are summarized by the single formula JTT 



SL _ I _ 1 (1 _ p.) 

2 6 3 V 3B/ /-s 

7 = 1 _ WE ' (5) 

1 3 B 

where P is pressure on the cold curve (T = equation of state), and B = —VdP/dV and 
B' = (dB/dV)/ (dP/dV) are the bulk modulus and its pressure derivative at T = 0. We 
will refer to Eq. (5) as the SDMVZ (Slater, Dugdale-MacDonald, Vashchenko-Zubarev) 
formula. Slater's derivation is based on the Debye-Griineisen definition (3) and assumes 
no volume dependence of the Poisson ratio. Dugdale and MacDonald used a simpli- 
fication of lattice dynamics in which the material is modeled as a lattice undergoing 
one-dimensional harmonic oscillations. Their derivation was improved by Vashchenko 
and Zubarev who considered three-dimensional oscillations of a lattice with interatomic 
interactions described by an anharmonic central potential. Equation (5) reduces to the 
Slater, Dugdale-MacDonald, and Vaschenko-Zubarev formulas for t — 0, 1, and 2, respec- 
tively. Eq. (5) has since been rederived by other researchers following different approaches 

rimimiE 
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An analytic model for 7(p) 



Melting curves and equations of state are usually based on simple functional forms for 
7(p). Dongquan and Wanxing |TJ|], for example, proposed the form •y(p) = joPo/p + 



2/3 (1 — po/p) s , which reduces to 7p ~ constant at small compressions and approaches 
2/3 as very large compressions; 7 does not differ from that given by Eq. (5) by more than 
10% for 1 < 5 < 3 fL5| . For 5 = 2 this form reduces to j(p) = 2/3 + 71/p + 72/ p 2 , a model 



frequently used for the construction of equations of state and melting curves || [lq| . 

At very high compressions, of order 10, a solid becomes a crystallized one-component 
plasma, i.e., a lattice of ions in a uniform neutralizing background of electrons [O. Several 



theoretical studies predict 7 = 1/2 for this limiting state of a solid. Kopyshev JT8[ 
calculated j(V) in the Thomas- Fermi approximation and found 7 — > 1/2 as V — > 0. 
Simple dimensional arguments by Hubbard |19[ also indicate that 7— >l/2asV A -^0. 



Additional theoretical studies that give 7 — > 1/2 include, but are not limited to, references 
PU| |2"T| , ^2] , p3[ . Nevertheless, some researchers assume 7— >-2/3asV— >-0 |24J because 



that is the limiting value of 7 as given by Eq. (5) with fixed t (see below). 

We now construct a simple, practical three-parameter model of 7(p) that accurately 
fits low-compression data, is equivalent to the experimental result that 7p' ?cff k, constant, 
Qcs > 1, for compressions up to 1.5, and limits to 1/2 as V — > 0. We assume that 
the Griineisen parameter is an analytic function of x = V 1 ^ 3 , essentially the interatomic 
distance, and that the coefficient of x in the Taylor- Maclaur in series expansion for 7 is 
non-zero. The simplest model that satisfies all of these requirements is 

1 (V) = ^ + c 1 V 1/3 + c 2 V q , ci, c 2 , q = const, q > 1, (6) 

Note that 7(V) — 1/2 is asymptotic to ciV 1 ^ 3 ; this is discussed in more detail in the next 
section. The term c 2 V q represents the contribution of the quadratic and higher-order 
terms in x which must sum to make 7 a concave-up function of V at small compressions. 
Eq. (6) can be rewritten in terms of density as 

1 7^ 72 

70) = o + + T> 7l ' 72 ' 9 = const > q>1 > W 

2 p 1 ' 6 p q 

Note that both 71 and 72 are dimensional parameters. 

Eq. (7) incorporates the low-compression power-law behavior 7p 9off w constant, ^eff > 
1, known from experiment. Experimental values of g e ff are typically in the range 1 to 2 p5|, 



but may be as high as 3 p6fl . The effective exponent can be defined by q e g = — d\wy /d\n p 
since 7 ~ p~ QcS ; in our model we have q c g = [7i/(3p 1 / 3 ) + q^/ P q ]/l(p)- Averaging over 
compressions from 1 to 1.25 we find ((fee) < 3 in general. Consider, for example, gold: 
with the corresponding parameters from Table 1, q e g varies from 4.7 at p = 19.3 g/cc to 
1.3 at p = 24 g/cc with an average value of 2.7. 

With 7(p) given by Eq. (7), the Lindemann equation (4) can be integrated to yield 
the melting curve 



T m (p) = T m (p r ) ^ exp 6 7l 




where p r and T m (p r ) are a reference density and corresponding melting temperature. 
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We now demonstrate that the asymptotic form of our model for the Griineisen param- 
eter follows from a generalization of the SDMVZ formula. 



Asymptotic 7(p) from a generalization of the SDMVZ 
formula 



Parshukov [E7| studied the compression dependence of 7 for lead, indium, and tin in the 



pressure range 0.6 to 4.0 GPa and found that it is almost identical for all three metals but 
is not described by Eq. (5) over the entire range of pressures for t a constant. He found 
that the Slater formula (t = 0) provided the most accurate fit at small compressions while 
the Dugdale-MacDonald one (t = 1) is optimal at higher compressions. He suggested that 
all three (t = 0, 1, 2) formulas be unified by making t compression dependent. Romain et 



al. |28[ noted that the melting curve of Al cannot be described in the Lindemann approach 
using the t — 0, 1, or 2 formulas for 7(p). Nagayama and Mori |^9| analyzed available 
experimental data on 16 metals and found that the data on 7 at P = are best described 
by the Slater formula, but at moderate compressions the data are best fit by the Dugdale- 
MacDonald formula which is approximately equivalent to 7/V 1 ' 1 = constant, close to the 
relation 7/V = constant that is routinely used in high-pressure studies. Nagayama and 
Mori noted that d In 7/d In V decreases with increasing compression from the value given 
by the Slater formula (t = 0) through that given by the Dugdale-MacDonald formula 
(t = 1) toward the value predicted by the Vashchenko-Zubarev relation (t = 2). An 
analysis of Hugoniot data shows that t follows a similar trend (0 — > 1 — > 2) in Al, Cu, 
and Ta These analyses of experimental data indicate that the parameter t is not a 



constant but a variable that increases with increasing compression. 

Irvine and Stacey |ff2| generalized the Vashchenko-Zubarev formula by accounting for 



non-central interatomic forces. They added a term / to the expression for the interatomic 
force constant and obtained the formula 

_ 2 6 _ 9B 18B 6 dP (n\ 
1 - 1 _ 4P J_ , l y J 

1 3B 3B 

which reduces to the Vashchenko-Zubarev formula for / = 0. By introducing the dimen- 
sionless parameter t via 

f = 2P(2-t), (10) 

Eq. (9) reduces to 

&_ _ 1 _ t f-t P_\ I P t y dt 

2 6 3 I 1 3B ) ^ 3B V dV /-, x 

7 = ; 2t p > ( n ) 

1 3 B 

which would coincide with Eq. (5) without the additional term (PV/3B) dt/dV in the 
numerator. Hence Eq. (11) is an extension of the SDMVZ formula, Eq. (5), to the case 
of density-dependent t. We note that Eq. (11) can be cast in the form 



1 b dlB-itP) 1 idlnlB-itP) 1 

™ — i L — i ± L (12) 

1 2 B -|tP dP 6 2 dlnV 6' K ' 
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which is equivalent to Eq. (5) if t is a constant. In view of the experimental evidence 
that t is a decreasing function of V we use Eq. (11) rather than (5) for our subsequent 
analysis. 

We consider Eq. (11) at ultrahigh pressures where the (T = 0) equation of state is 
accurately given by |32| 



P = aV- 5/3 e~ bVl/3 , a, b = const > 0. (13) 
Eq. (13) includes an exponential screening correction to the equation of state of a free 

O 

electron gas, namely P = aV~ b ' 3 where a = 2.337 Z 5 ' 3 TPa A 5 , Z being the atomic 
number JJIJ. It agrees with very accurate numerical Thomas- Fermi-Dirac results to within 
2% over the compression range of 1 to 15 with the exception of the alkali and alkaline-earth 



metals |3TJ. Using the equation of state (13) in Eq. (11) we obtain 



_ 4(5 - 2 t) + 2(4 - t) b V 1 ' 3 + b 2 V 2 / 3 + 6V dt/dV 
1 ~ 6(5 - 2 t + & V 1 ' 3 ) ' ( ' 

It follows that 7-^2/3 as V — ► for every asymptotic value of t except t = 5/2. If t — 5/2 
then 7 — > 1/2 as V — > 0, in agreement with theoretical predictions. We conclude that t is 
always asymptotic to 5/2. Recent computer calculations of the compression dependence 
of t for 7-Fe [13] corroborate this conclusion and also show that t can saturate at quite 
moderate compressions: t = 2.46 ± 0.03 at compressions of only 1.5 to 2. 

The cold-curve (T = 0) pressure (13) is an analytic function of x everywhere except 
at the origin (infinite compression) where it has a pole of order five. (The majority of 



model equations of state are analytic in x except for poles at x — 0. See Holzapfel ]32 
for a list of 16 examples.) Consequently, both P/B and B' are analytic for x > 0. Given 
the analyticity of 7, P/B, and B' for x > 0, it follows from Eq. (11) that t must also be 
an analytic function of x for x > 0, hence it can be represented by the power series 

00 

t = b/2~Y J Ux i . (15) 

i=l 

The first sub- leading coefficient, tj, must be non-negative if t is in fact a monotonically 
increasing function of compression. Substituting this series in Eq. (14) and expanding we 
find that the asymptotic form of 7 is 

2 60 + 12 ti 

It is expected that the V 1 ! 3 term is always present in the series expansion of 7, i.e., 
b 2 + 2 b t\ — 2 1 2 7^ 0. Thus the generalized SDMVZ formula leads to 7 ~ 1/2 + cix, x — > 0, 
in agreement with Eq. (6). 

In summary, we have (i) generalized the SDMVZ formula for 7(p) to account for 
a density- dependent t, (ii) shown that 7 goes to its theoretical high-pressure limit 1/2 
only if t is asymptotic to 5/2, and (iii) demonstrated that our model for the Griineisen 
parameter and the generalized SDMVZ formula, Eq. (11), have the same asymptotic 
behavior provided t is an analytic function of x, and t(0) = 5/2. 
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Determination of model parameter values 



In this section we outline the procedure for extracting the numerical values of the three 
model parameters 71, 72, and q from data and determine their values for 20 elemental 
solids. 

Ideally, one would fit our model (7) to data on 'y(p) and extract the values of 71, 72, and 
q, but such data sets are very rarely available. Alternatively, one could fit the functional 
form (8) to T m (p) data, but such data are very rarely available either since experiment 
usually determines T m (P). Typically, the only available data on 7 are 7(^300 = p{P = 
0, T = 300 °K)), which is approximately equal to jth(T = 300 °K), and 7(p m = p(P = 
0, T = T m )) which can be found by equating dT m /dp at P = as given by the Kraut- 
Kennedy law |33f 



Tm \P) T m (p r 



1 + 2 Up 




;i7) 



to dT m /dp = (B m / p m ) dT m /dP, where B m and dT m /dP at P = are taken from ex- 
periment. The value of B m is obtained by calculating bulk moduli at low temperatures 
from measured single- crystal elastic constants and then extrapolating to T m . The pres- 
sure derivative dT m /dP is obtained either directly from low-pressure melting curve data 
or indirectly from isobaric-heating measurements of AH (H is the enthalpy) and AV 
across the melting transition so that dT m /dP = T m AV / AH (Clausius-Clapeyron equa- 
tion). Eq. (17) is obtained by expanding T m (p) in a power series about p m using Eq. (4). 
Equation (7) with p = p 300 and p = p m provides two conditions needed to determine the 
three parameters 71, 72, and q. The third condition comes from the ultrahigh pressure 
limit, which we discuss next. 

The melting curve of a solid at ultrahigh pressures is described by the equation 

Z 2 e 2 

T m , (18) 



akBT n 



where a = {3v/Att) 1 ^ is the Wigner-Seitz radius (v being the Wigner-Seitz volume) and 
T m , a dimensionless constant, is the OCP coupling parameter at melt. The value of 
the coupling parameter is 170 — 180 for a body-centered cubic (bcc) OCP crystal |34| 
(the recent calculation of T m for a bcc crystal by Potekhin and Chabrier []35] gave T m = 
175.0 ± 0.4), and as high as 200 - 210 for a face-centered cubic (fee) OCP crystal |fj, |7|. 
In the following analysis we make no distinction between bcc and fee OCP crystals and 
take T m = 180. 

It follows from Eq. (18), and Eq. (8) in the limit p —>■ 00 with p r = p m that 

,i/3 „J^M = r 4 ^ 1/s e2 22 



r m wr«p + =It rr- (19) 

I p m nP-m ) \ o / Kb 1 r 



in 



where v m , the Wigner-Seitz volume at the zero-pressure melting point, equals a m /2 for a 
bcc and a^/4 for a fee crystal; a m is the lattice constant. 

We obtain the values of the parameters 71, 72 and q by simultaneous solution of three 
non-linear equations, namely Eq. (7) for p = p 300 and p = p m , and Eq. (19) with T m = 180. 
The value of 7(^300) , which is very rarely measured directly, can be approximated by 
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7</i(300) = aBx/ipCv) where a, Bt, and Cy are all measured at 300 °K; the necessary 



room-temperature data can be found in ref. [38], and the data on 7^(300) in ref. [39 
Parameter values for 20 elements are shown in Table 1; 71 and 72 are given to three 
significant figures, and q to two significant figures. 

The uncertainty in the value of 7(p m ) is determined by the error bars on the measured 
melting temperatures and densities, usually a few percent (except perhaps for laser-heated 
diamond anvil cell (DAC) melt data, see below). Typical uncertainties in a, Bt, and Cy 
result in an uncertainty of roughly 5% in the value of 7(p3oo) ( we consider the difference 
between 7(^300) and 7^(300) to be negligible). These uncertainties in 7(p m ) and 7(^300) 
imply an uncertainty ~ 5% in 71, values of q that are accurate to ~ 40%, and order-of- 
magnitude uncertainties in 72 that are attributable to a very strong sensitivity to the value 
of q. Let us take, as an example, copper; with y(p m ) — 2.48 and 7^(300) = 2.19 ±0.1 we 
find, neglecting uncertainties in the density, 71 = 1.84 ± 0.08, 72 = 430 — 1.95 • 10 6 , and 
q = 4.7 ±2.0. 

Although the values of 72 have order-of-magnitude uncertainties, we note that if Eq. 
(7) is rewritten in a "scaled" form, 7 = 1/2 ± 7i(p3oo/p) 1//3 + 72(p3oo/p) 9 , the values of 
both 71 and 72 are of order 1, as can be seen in Table 1 where we also include the values 
of both 71 and 72 to three significant figures. 



element 


Z 


7!, 


72, ( 


JL\<? 

cc 1 


71 


72 


q 


P300, g/cc 


Pm, g/cc 


T m (p m ), °K 


Ne (fee) 


10 


1.04 


5.17 




0.91 


2.10 


2.2 


1.507 


1.435 


24.6 


Mg (hep) 


12 


0.79 


4.33 




0.66 


0.53 


3.8 


1.740 


1.640 


923 


Al (fee) 


13 


0.84 


45.4 




0.60 


1.40 


3.5 


2.700 


2.550 


933.5 


Ar (fee) 


18 


1.19 


16.3 




0.90 


1.86 


3.8 


1.771 


1.622 


83.8 


Fe (bee) 


26 


1.72 


2.66 


■ 10 3 


0.86 


0.25 


4.5 


7.870 


7.270 


1811 


Co (fee) 


27 


1.81 


6.28 


• 10 4 


0.88 


0.39 


5.5 


8.830 


8.180 


1768 


Ni (fee) 


28 


1.85 


5.60 


• 10 5 


0.89 


0.38 


6.5 


8.900 


8.220 


1728 


Cu (fee) 


29 


1.87 


2.31 


• 10 4 


0.90 


0.78 


4.7 


8.930 


8.370 


1357.7 


Zn (hep) 


30 


1.91 


1.84 


• 10 3 


0.99 


1.05 


3.8 


7.140 


6.900 


692.7 


Mo (bee) 


42 


2.06 


1.40 


• 10 6 


0.95 


0.19 


6.8 


10.21 


9.650 


2896 


Pd (fee) 


46 


2.40 


3.34 


• 10 6 


1.05 


0.25 


6.6 


12.02 


11.29 


1828 


Ag (fee) 


47 


2.23 


9.63 


• 10 4 


1.02 


1.21 


4.8 


10.49 


9.850 


1235.1 


Cd (hep) 


48 


2.43 


2.47 


• 10 4 


1.18 


0.97 


4.7 


8.650 


8.420 


594.2 


In (bet) 


49 


2.43 


6.14 


• 10 3 


1.25 


0.80 


4.5 


7.310 


7.200 


429.8 


Sn (bet) 


50 


2.37 


2.37 


• 10 3 


1.22 


0.83 


4.0 


7.300 


7.200 


505.1 


Pt (fee) 


78 


3.21 


1.13 


•10 11 


1.16 


1.01 


8.3 


21.45 


20.19 


2041 


Au (fee) 


79 


3.21 


1.97 


•10 12 


1.20 


1.62 


9.4 


19.30 


18.29 


1337.6 


Tl (bee) 


81 


3.17 


3.74 


•10 7 


1.39 


0.54 


7.3 


11.85 


11.55 


577 


Pb (fee) 


82 


3.09 


8.21 


• 10 8 


1.38 


0.89 


8.5 


11.34 


11.05 


600.6 


U (bee) 


92 


3.39 


1.05 


•10 11 


1.27 


0.57 


8.8 


19.05 


17.60 


1405 



Table 1. Numerical values of the parameters entering Eqs. (7) and (8) for 20 elemental 
solids. The crystal structure indicated for an element is that from which it melts at zero 
pressure. For Ne and Ar the entry under p 30 o is the value of p(T = 0). 



Melting curves 

In this section we calculate melting curves using Eq. (8) with parameters determined in 
the previous section and compare the results with the available experimental data. 
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In Fig. 1 we compare our theoretical aluminum melting curve to experimental data ||0] . 
Although the theoretical melting curve looks like a fit to the datapoints, the parameters 
7x, 72, and q were not obtained from such a fit but rather from the zero-pressure data 
7 (p 3 oo = 2.7 g/cc)=2.5 H|, 7 (p m = 2.55 g/cc)=2.83 H, and T m (p m ) = 933.5 °K, and 



.5 m, 7(p m = 2.55 g/cc)=2.83 
Eqs. (7) and (19). In Figs. 2-4 we compare our theoretical melting curves for argon, nickel, 
palladium, and platinum to melting data available in the literature. Agreement between 
our theoretical curves and experiment is good with the exception of the lower-pressure 
data on nickel which were obtained from a laser- heated DAC f41 |. This discrepancy may 
be due to systematic errors in the DAC data 



The Z dependence of 71 

It is evident from Table 1 that 71 is a slowly increasing function of atomic number. We 
have fit the forms 7 x = C\ • Z n + C2 and 71 = G\ • (In Z) n + C2, where C±, C2 and n are 
constants, to the table entries and find that the minimum x 2 , 0.11, is realized for both 
71 = (12/11) Z 1 / 3 - 11/7 and ^ = (2/21) (InZ) 7 / 3 + 2/13. These fits may be used to 
predict 71 in those cases where data provide only two constraints on 71, 72 and q. In Fig. 
5 we plot 71 = (2/21) (InZ) 7 / 3 + 2/13 along with the 7i(Z) entries in Table 1. 

We note that the 71 entries in Table 1 can be fitted to the same functional form, i.e., 
71 (Z) = (7/26) Z 1 / 3 + 7/59, but the accuracy of the fit is much lower than that of the 
7i(Z) fit: x 2 = 0-97. 



Concluding remarks 

Our model for 7 accurately fits low-pressure data and agrees with theoretical predictions 
that 7 — > 1/2 at ultrahigh pressures. Its accuracy cannot be determined at intermediate 
compressions because there are no experimental data. However, comparison of our calcu- 
lated melting curves for aluminum and argon to the corresponding data in Figs. 1 and 2 
demonstrates that our model is accurate up to compressions of at least 2. (Our calculated 
melting curve for copper, not discussed in this paper, is in excellent agreement with a new 
SESAME copper melting curve || up to compressions ~ 10.) In our model the deviation 
of 7p 9cft from a constant over the range of compressions from 1 to 1.5 is generally less 
than 10%. Finally, we note that our model for 7 can be used to calculate shear moduli 
along the solidus using Eq. (8) and the formula G(p,T m (p))/(p ■ T m (p)) = constant p3|] . 
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Figure captions 



Fig. 1. Melting curve of Al: Eq. (8) 



with the Al parameters from Table 1 vs. the exper- 

The experimental 



imental data |40| and the results of calculations (larger points 
points have error bars of ~ 100 — 150 °K which are not shown on the plot. 



Fig. 2. Melting curve of Ar: Eq. (8) with the Ar parameters from Table 1 vs. the experi- 
mental data |45| and the results of calculations jSJ (larger points). 



Fig. 3. Melting curve of Ni: Eq. (8) with the Ni parameters from Table 1 vs. the ex- 



perimental data [41 . The upper three points are shock-melting data from ref. [47 which 



does not quote error bars. 



Fig. 4. Melting curves of Pd and Pt: Eq. (8) with the Pd and Pt parameters from 
Table 1 vs. shock-melting data ||8| (upper two points of each curve, no error bars given) 
and zero-pressure datapoints (lowest points on each curve). 

Fig. 5. Comparison of 71 = (2/21) (lnZ) 7 / 3 + 2/13 to the 7i(Z) entries in Table 1. 
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